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Abstract 

We propose a new method to evaluate the chiral charge flux, which is converted into baryon 
number in the framework of the charge transport scenario of electroweak baryogenesis. 
By the new method, one can calculate the flux in the background of any type of bubble 
wall with any desired accuracy. 
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1 Introduction 



In the scenario of electroweak baryogenesisjj]], the chiral charge flux, which is converted 
to baryon number by the sphaleron process, is one of key ingredients. That is supposed to 
be generated by two mechanisms, both of which depend on the interaction of the fermions 
with the CP-violating background of the bubble wall created at the electroweak phase 
transition (EWPT). One of the mechanisms is the charge transport scenario, which is 
quantum mechanical and effective for thin walls, the other is the spontaneous baryogenesis, 
which is classical and effective for thick walls. 

As for the charge transport scenario, the flux is determined by the difference in the 
reflection coefficients of the left- and right-handed fermions, as well as by that in the 
distribution functions in the broken and the symmetric phases divided by the bubble wall. 
The reflection coefficients are calculated by solving the Dirac equation in the background 
of the bubble wall, which is accompanied spatially varying CP violation. Cohen, et a/.|2|] 
first evaluated the flux by solving the Dirac equation numerically, assuming the kink- 
type profile for the modulus of the Higgs field and a function linear in the kink for the 
phase. On the other hand, we formulated a perturbative method applicable for small 
CP-violating phase and proved various relations among the reflection and transmission 
coefficients ||, |J. Although it enables us to derive analytical relations, it can actually be 
applied only to the case in which the unperturbed Dirac equation is exactly solved and 
the CP-violating phase is sufficiently small. However the profile of the bubble wall should 
be determined by the dynamics of the EWPT. We pointed out, by solving the equations 
of motion for the gauge-Higgs system, that the CP- violating phase can be of 0(1) even if 
that is small in the broken phase at T = 0, according to the parameters in the effective 
potential ||. For a larger CP- violating phase, we naively expect larger baryon number to 
be generated. So we need a method to evaluate the chiral charge flux for arbitrary profile 
of the bubble wall. 

In 0, the one- dimensional Dirac equation is solved within a finite range [0, Zq] includ- 
ing the bubble wall, imposing the boundary condition in which plane waves are injected 
at the boundary. Strictly, a plane wave is not an eigenfunction at any finite z, so that 
the choice of the range will affect the precision of the results. In this paper, we pro- 
pose another procedure to solve the Dirac equation and to calculate the reflection and 
transmission coefficients. Although the method is numerical, its precision can be so con- 
trolled that one can obtain the results with any accuracy. In section 2, we write down the 
CP-violating Dirac equation and give the expressions of the reflection and transmission 
coefficients of the chiral fermions. The procedure of the numerical analysis is given in 
section 3, and it is applied to some profiles, including those which we found in || and 
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that in [^] for comparison, in section 4. The final section is devoted to summary and 
discussions. 

2 CP- Violating Dirac Equation 



2.1 The Dirac equation 

We consider a fermion in the background of a complex scalar, to which it couples by a 
Yukawa coupling. For the scalar to have a nontrivial complex vacuum expectation value 
( VEV) , we assume it is a part of an extended Higgs sector of the electroweak theory, such 
as MSSM or more general two-Higgs-doublet model. At the electroweak phase transition, 
the Higgs changes its VEV from zero to nonzero value near the bubble wall created if the 
phase transition is first order. 

Then the Dirac equation describing such a fermion is 

i<fh\){x) — m(x)P R i/;(x) — m* '(x)PliI>(x) = 0, (2.1) 

where P R = -4pk,Pi, = and m(x) = —f((p(x)) is a complex-valued function of 

spacetime. Here / is the Yukawa coupling. When the radius of the bubble is macroscopic 
and the bubble is static or moving with a constant velocity, we can regard m(x) as a static 
function of only one spatial coordinate: 

m{x) = m(t, x) = m(z). 

Putting 

^(x) = e M-£*+PT-x T ) ; ( a = +loT _ X ) (2.2) 

Q2.1|) is reduced to 

a(7°£ - 7tPt) +h 3 d z - m R {z) - i^m^z) iPe{~Pt, z) = 0, (2.3) 

where 

Pt = (p\p 2 ), x T = (x 1 ,x 2 ), Pr=|pr|, ItPt = 7 V + tV> 
m(z) = m R (z) + imi(z). 

If we denote E = E* cosh 77 and px = E* sinhr] with E* = E 2 — p T , Pt in (|2.3| ) can be 
eliminated by the Lorentz transformation x/j 1— > ip' — e~ r?7 ° 75 '?/': 

d z ^ E (z) = i7 3 [-<x£*7° + m R (z) + i^m^z)] ^j e (z). (2.4) 



3 



In our work on the perturbative treatment of CP violation, we used the eigenspinor of 7 3 , 
since unperturbative part of the Dirac equation contains only 7 s . Here we shall use the 
chiral representation of the 7-matrices as in 0. So the expression of the Dirac equation 
is the same as that in 0, but we write it for self-containedness. 
In the chiral basis, the 7-matrices are represented as 
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-a k 
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If we write the four-spinor as 



the Dirac equation is decomposed into two two-component equations: 



d 



ip 4 (z) 



aE* m*{z)\ fipi(z) 

-m(z) -aE*)\^(z) 

aE* 777(2) \ (^a{z 

-m*{z) -aE* ) \fa(z 



(2.5) 
(2.6) 



Now let us introduce dimensionless parameters. Suppose that a is the parameter 
of mass dimension, whose inverse characterizes the width of the bubble wall. We shall 
comment on practical choice of it in the later section. Then put 



x = az, 
t = E*/a, 



(2.7) 



i(x)e ie{x) = m{z)/a 

Since we are interested in the bubble wall background, the m(z) is supposed to behave as 

771(2) — > m e ie °, as z — >• +00, 
777,(2) — > 0, as z — > —00, 

where m is the fermion mass and 8q is the CP-violating angle in the broken phase at the 
phase transition. We also put 

£ = m /a. (2.8) 
In terms of these dimensionless parameters, the Dirac equation is expressed as 



Or 



0,- 



V>i(x) 

4> 3 (x) 

ip 4 (x) 



at 








£(x)e ie 


W -at ) 
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m _ ae ) 
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(2.9) 
(2.10) 



4 



This is the equation we shall numerically analyze. 

Before investigating the Dirac equation, let us write down the chiral currents in our 
basis. They are defined by 

f L {x) = ^(x)yl^^(x), 

f R (x) = ^ x )^ l -±pL^ x ). 

It is easy to see that the relevant currents are written as 

f L (x) = ^(aOI 2 -!^)! 2 , (2.11) 
f R (x) = |^(x)| 2 -|^(x)| 2 . (2.12) 

The asymptotic forms of these currents will be used to calculate the transmission and 
reflection coefficients. 



2.2 Extraction of the oscillating factor 

Here, for definiteness, we consider a scattering state in which the incident wave is coming 
from the symmetric phase region (x ~ — oo) so that there is only the transmitted wave 
(right moving) in the broken phase (x ~ +oo). This boundary condition is expressed as 

ip a (x) oc e aax as x — > +oo, 

where a = i\/e' 2 — £ 2 and a denote the spinor index. Our aim is to find out the coefficients 
of the right and left movers in the symmetric phase region of this state, from which we 
can evaluate the chiral transmission and reflection coefficients. 

For this purpose, we numerically integrate the Dirac equation starting from the plane 
wave in the deep broken phase, x = +oo. However it is rather difficult to implement this 
situation numerically, since preparing the plane wave at x = oo is impossible, while it 
is the eigenstate only there. One may think that an appropriate change of the variable 
which brings x = +oo to some finite coordinate makes the wave function rapidly oscillate 
near the boundary of the region. To remedy this, we factorize the oscillating part of the 
wave function. 

Before performing the factorization, let us remove the effect of CP violation in the 
broken phase by a unitary transformation of the two-spinor. Introducing a unitary matrix 




(2.13) 



and transforming the spinor as 
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the Dirac equation is reduced to 

R (4>i(x)\ f (f(x) -i[ae + £g(x)]\ (i>i(x)\ . , 

x \Mx)J {-ilae-Zgix)} -£/(*) JU(i)J' 1 ^ 

a (M x )\ ( -£f( x ) -i[vt + £g{x)] \(Mx)\ . . 

x \Mx)) ^ \-i[ae-ig{x)] £f(x) ) \fa(x) J ' ^ iDJ 

where 

Cf(x) = £(x) cos(#(x) - 6 ), ig{x) = £{x) sin(0(x) - 6 ). (2.17) 
These satisfy the boundary conditions 

f(x) — > 1, g(x) — * 0, as x — ► +oo 
f(x) — > 0, g(x) — >■ 0, as x — * — oo. 

In the following, we consider only the positive frequency wave (er = +) for definiteness. 
For the boundary condition we concern, the wave functions behave as, at x ~ +oo, 

Mx) ~ C a e ax . (2.18) 

The asymptotic form of ( |2.15| ) yields 

(f - a)de ax = ieC 3 e ax , (£ + a)C 3 e QX = -ieC ie QX , 

which gives 

C 3 = ^d. (2.19) 

Similarly we have 

C 4 = ^^C 2 . (2.20) 

Now we factorize the wave function as 

Mx) = C 1 e ax Xi(x), 
Mx) = C 2 e ax X 2(x), 

Mx) = L^de-XsCar), (2 ' 21) 

Mx) = ^-^C 2 e ax X i(x). 
te 

The boundary condition is now 

Xa {x = +oo) = 1. (2.22) 
In terms of Xa, the Dirac equation (|2.15|) and (|2.16|) is expressed as 



x'i(x) = (C - a)(xi(a;) - Xs(^)) + f (/(a?) - l)xi(a?) g(x)xs(x) 

Xz{x) = (£ + a)(Xi(aO - Xs(a;)) - £(/(a0 - l)Xs(^) - ^j-^- <y(z)xi(g) 
XaO*) = (6 - «)(X2(x) - X4(x)) + f (/(ar) - l)x 2 (x) + ^ a ^ (x)x 4 (^) 

xK^) = (f + «)(x2(^) - X4(^)) - f (/(z) - i)x4(») + ^^ + a ^ ( x ) X2 ( x ) 



(2.23) 



(2.24) 



Note that the right-hand sides of the above equations vanish at x = +00. 



2.3 Reflection and transmission coefficients 



To determine the transmission and reflection coefficients, we need their asymptotic forms. 
Suppose that the wave function behaves as, at x ~ —00, 



where f3 = ie. Then the incident and reflection currents are given by 



(2.25) 



refl 



W-\A Z \\ 



refl 



(2.26) 



\B A \*-\B 3 \<, = \B^-\B 2 \\ 

In the asymptotic region (x ~ —00), the coefficients A a and _B a are given by 

e -Px / 1 \ 

= —ha(x) + -^ a (x)\, 

e px / ! \ 

#a = -y ( 1>a(x) ~ p1p a (x) J . 

Strictly speaking, these are not independent of x at any finite x. However, we shall use 
these expressions to calculate the transmission and reflection coefficients for so large x 
that they stay almost constant. 



Because of (|2.14|) 



^4 
"02 
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4>3 
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( 7^2-^4) 
1^2 + ^4 



In terms of Xa{x), the "constants" are expressed as 
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\A 4 



\Ci\ 



ICil 
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I Co 



1 

(l + -g)Xi(aO + _ 

a 1 
(l + -)Xi(x) + -x'i(^) + 

a 1 

(1 + g)X2(x) + gX^) + 

(1 + ^)X2(X) + -X2(^)- 
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/3 
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P 




— a 



1.7 



(1 + ^)Xa(x) + ^XsO*) 

a 1 
( 1 + ^)X3(x) + -x' 3 (x) 

a 1 

(1 + -g)X4(a:) + -x*(x) 

a 1 
(1 + -^)Xa(x) + -X 4 (x) 



(2.27) 
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and 



\Bi\ 
I-Ssl 

\Ba\ 
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(1- 


j)Xi{x) 


- \x^) 
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%)Xa{x) 


1 

'p 



-Mx) 



- 13Xii.x) 



(2.28) 



where ' denotes the derivative with respect to x. On the other hand, the transmitted 
currents are 



(^) 



trans 



trans 1 , ,. 

= 5 IC|- 



£ — a 



P 
£ — a 



1 + 
1 + 




£ - a 







(2.29) 



Equipped with these coefficients, we are ready to evaluate the transmission and reflection 
coefficients. 

The transmission and reflection coefficients of the chiral fermions are defined by 



L^L(R) 



Ohm) 



T, 



R— >L(_R) 



(J'l) 1 



trans 



R 



R(L)^L{R) 



(jl) mC 



( \ inc 
i£) =0 



refl 



3r(l) 



(2.30) 



(2.31) 



(2.32) 



(V V 



The condition 



Jl(r) 



is realized by appropriately choosing the ratio of C\ and C 2 . 



For example, to realize {j%) mc = 0, we take 
ICil 2 



en I £ — (x 



<i + f)xi(*) + ±xl<*)-^ 



a 1 
(1 + ^)x 4 (*) + -^(z) 

a 1 

(1 + -g)X3(aO + -pXsix) 



(2.33) 



so that (j 



|Ai| 2 - |A 2 | 2 = 0. Given these |C;| 2 , (gj^ ~( ggD yield the chiral 



transmission and reflection coefficients by use of (|2.26| ) and Q2.29] ). 
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3 Numerical Analysis 

To realize the boundary condition ( |2.22| ) at x = oo on a computer, we change the variable 
of infinite range to another of half-infinite range. For definiteness, we take the following 
new variable 

y = e~ az = e~ x G (0, oo). (3.1) 
This transforms ^ i— > — y-^ so that it brings the singularity at y — into the Dirac 



equation. As we show below, the right-hand sides of ( [2.23|) and ( [2.24|) well behave to 
cancel such singularity. 

The bubble wall connecting the broken and the symmetric phase is expected to behave 

as 

p(z)~l + e" b2 , as z -> +oo , , 

f \ b'z V-*-^) 

p(z) ~ e , as z — > — oo 

where 6(> 0) and b'(> 0) will be order of the Higgs masses. Hence if we choose a to be 
smaller than the Higgs mass (a < b), the potentials in the Dirac equation behave as, at 
y ~ (x ~ +oo), 

+ g(y)~y u , 

with [M, v > 1. Now put 

oo oo 

f(y) = l + E fnv\ g(y) = y u Y. 9 n y n . (/„ ^ o, g + o) 

n=0 n=0 

To find the asymptotic behaviors of Xi(y) an d X^{v)i expand them as 



Xi (y) = l+y p J2 a nV n , (R 3 P > 0, a ^ 0) 



n=0 

oo 



* 3 (z/) = i + y a E & ^ n - (R 3 ^ > o, 6 7^ o) 

n=0 



(3.3) 



Then (2.23) in terms of y yield 



-EG°+ «Kf n = (£ - «) E a ^ n - E + £ ^ p E fny n (i + y p E a n?/ n ) 

n=0 \n=0 n=0 / n=0 n=0 

-WkZ^ltf-P £ ^(i + ^ £ ^v"), (3.4) 

e n=0 n=0 

oo / oo oo \ oo oo 

-J2(v + n)b n y n = ^ + a)(y^Y,^y n -T, b ny n )-^y fl " 7 Efny n a + y a E b ny n ) 

n=0 V n=0 n=0 / n=0 n=0 

^V"* E ^ n (l + 2/ P E (3.5) 

6 n=0 n=0 
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When p < p, a or v < p, a, the right-hand sides of the both equations have singularities. [] 
When p > p,a and v > p,a, we must have p = a. This is because, otherwise, the right- 
hand side of either (|3.4j) or (|3.5|) has a singularity, boy a ~ p or a y p ~ <T , which cannot vanish 
because of a ^ and b ^ 0. Then 0(y°) terms give 

(£ - a + p)a - (£ - «)&o = 0, 
(£ + a)a - (£ + oi — p)b = 0. 

To have nontrivial ao and bo, 

(£-a + p)(Z + a-p)-(Z-a)(Z + a) = p{2a - p) = 0, 

which has no real solution p > 0. Therefore to avoid any singularity in the right-hand 
sides, we are left with two possibilities: 

p — a — p < v, or p = a = v < p. (3.6) 

In these cases, the lowest order equations determine a and b in terms of f and/or go- 
Needless to say, the same argument also applies to the remaining components of the Dirac 
equation. 

Thus the Dirac equations in terms of y are free from any singularity, so that they 
can be numerically integrated by use of a standard algorithm such as the Runge-Kutta 
method. Now we summarize our procedure to calculate the chiral reflection coefficients: 

1. Prepare the potential f(y) and g(y)- These would be determined by the model we 
choose through the equations of motion. 

2. Integrate the Dirac equation 

-v dj ^- = (£ - «)(xi(y) - Mv)) + t(f(v) ~ l )xM - ^^(v)xs(v), 
-y^l = (£ + a)(xi(y) - xM) - C(f(y) - i)xs(y) - ^ ± ^-g(y)xi(y), 

dy e 

_ y ^M = £ + a)iX2iy) _ Uy)) _ afiy) _ l)Uy) + i(i±^)^ (y)X2(y)) ° 

ay e 
starting from the boundary condition Xa(y — 0) — 1 to large ?/ with a step Ay. 
Within each step, one can control the precision by choosing sufficiently small meshes 
or using quality control algorithm which automatically divides the mesh until the 
desired precision is reached. 



1 One may think that in ( fc.4|) when y — v < p, they cancel with each other if e/o = (£ — a) go- But 
this is impossible for real fo and go- 
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3. At the end of each step, evaluate the reflection and transmission coefficients. 



4. Proceed to larger y, until the relative error in the coefficients comes to be smaller 
than the prescribed precision. 

We shall see how this program works in the next section. 

4 Examples 

Here we take the following potentials to illustrate the calculation of the chiral reflection 
coefficients. The chiral charge flux can be obtained by the procedure outlined in |4[] 
integrating the reflection coefficients multiplied by the statistical distribution functions of 
the fermions. 

, 1 — tanhx 



[i) e(x)e^)=e 1 + t o anh3: exp 



iAO- 



This with A9 = —tt was taken in 0. 

(h) e(x)e^) = e 1 + t 2 anhx e^), 

with 9{x) being the solution to the equations of motion of the gauge-Higgs system, 
depicted in Fig. 2 of ||, satisfying 9 = 9{y = 0) = 1 and 9(y = oo) = 0. 

(hi) e(x)e^)=e 1 + t 2 anh3: e^), 

with 9{x) being the solution to the equations of motion of the gauge-Higgs system, 
depicted in Fig. 3 of ||, satisfying 9(y = 0) = 9{y = oo) = 0. 

In the charge transport scenario, the generated baryon number density is given by || 

s ~ " ttV Kaw uT* ' (4 - lj 

where M is the model-dependent factor of 0(1), g* is the radiation degrees of freedom 
at T c and is of O(100), notyy is the sphaleron transition rate in the symmetric phase and 
is of O(10~ 6 ) and u is the wall velocity, r stands for the transport time, which is given 
by the average time the reflected fermion spends in the plasma before caught up with 
the bubble wall. It is expected to be of O(10) ~ 0(1000) times the thermal correlation 
lengths, depending on the wall velocity, if the forward scattering is enhanced. In the 
original scenario, the weak hypercharge was taken as the chiral charge injected into the 
symmetric phase. In the following, we mean by the chiral charge flux the dimensionless 
quantity normalized as 

Fq (4.2) 



uTZ{Q L -Q R 
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case (i) 

In terms of the variable y, the potentials are written as 

Now we take A9 = — tt. To show how the program in the previous section works, we 
calculated the reflection and transmission coefficients of the left-handed fermion for £ = 1 
and e = 1.2. The result is listed in Table [I] for first several y. The column 'unitarity' 



y 


Rl^r 




T L _ L + T L . 






unitarity 


20 


5 


42508153 


X 


10- 


-1 


4 


57493406 


X 


10" 


-1 


1 


00000156 


40 


5 


42731549 


X 


10- 


-1 


4 


57268549 


X 


10- 


-1 


1 


00000010 


60 


5 


A21 


319216 


X 


10- 


-1 


4 


57180804 


X 


10- 


-1 


1 


00000002 


80 


5 


42^ 


330565 


X 


10- 


-1 


4 


57169442 


X 


10- 


-1 


1 


00000001 


100 


5 


42^ 


326265 


X 


10- 


-1 


4 


57173738 


X 


10- 


-1 


1 


00000000 


120 


5 


A21 


320329 


X 


10- 


-1 


4 


57179672 


X 


10- 


-1 


1 


00000000 


140 


5 


A2i 


315576 


X 


10- 


-1 


4 


57184425 


X 


10- 


-1 


1 


00000000 


160 


5 


42£ 


312238 


X 


10" 


-1 


4 


57187764 


X 


10- 


-1 


1 


00000000 


180 


5 


42^ 


310021 


X 


10- 


-1 


4 


57189980 


X 


10- 


-1 


1 


00000000 


200 


5 


42^ 


308603 


X 


10- 


-1 


4 


57191398 


X 


10- 


-1 


1 


00000000 


220 


5 


42^ 


307730 


X 


10- 


-1 


4 


57192270 


X 


10- 


-1 


1 


00000000 


240 


5 


42^ 


307222 


X 


10- 


-1 


4 


57192778 


X 


10- 


-1 


1 


00000000 


260 


5 


42^ 


306955 


X 


10- 


-1 


4 


57193046 


X 


10- 


-1 


1 


00000000 


280 


5 


A2t 


306844 


X 


10- 


-1 


4 


57193157 


X 


10- 


-1 


1 


00000000 


300 


5 


421 


306833 


X 


10" 


-1 


4 


57193167 


X 


10" 


-1 


1 


00000000 



Table 1: The reflection and transmission coefficients of the left-handed fermion. These 
should become constants for y ~ 00. 

means the sum of the two columns left to it. This shows that to get the precision of 
four digits, it is sufficient to integrate the equations up to y = 100, which is amount to 
x = — In 100 ~ —4.6052 and is regarded far from the bubble wall. In this case, one can 
obtain the same data of AR = Rr-^l — Rl^r as a function of (e, £) as shown in Fig. 2 of 

§• 

For A# = 0.001, we confirmed that our numerical method yields the same AR as that 
obtained by the perturbative method in ||. For A# = —1, which is equivalent to the 
profile 

(4.4) 



e(x)e^)=£ 1 + t 2 anhX exp 



1 + tanh x 
i 



up to the unitary transformation ( |2.14| ), we calculated AR by the numerical method 
developed here. The chiral charge fluxes for various mass and wall width are depicted in 
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Fig. 1. 
case (ii) 

In this case, 6q = 1 and the maximum of g(x) is about 0.4, which is comparable with f(x) 
at the same x. So we did not evaluated the chiral charge flux in |J, using our perturbative 
method[|J. This profile, which is a solution to the equations of motion to the gauge-Higgs 
system with the kink-type p(x), indicates that 9{x) is no longer proportional to the kink 
when #o is 0(1). This should be compared with the profile (|4.4j) , which satisfies the same 
boundary conditions. We performed the numerical calculation of the chiral charge flux, 
which is given in Fig. 2. This shows that this profile slightly enhances the flux over the 
profile employed in [0 by about a factor of 2. The former is about 80% larger near the 
peak, and decreases more slowly for large mass or thick wall than the latter, 
case (iii) 

Although the maximum of 6(x) reaches 0.37, the potential f(x) has a kink-like shape and 
g(x) is sufficiently small as shown in Fig. 4 of 0. We expected that the perturbative 
method is applicable to this case and evaluated the chiral charge flux, depicted in Fig. 5 
of ||. To confirm this, we calculated AR for several (e, £). One of the results is given 
in Fig. 3 for £ = 0.4 and m = 1. This shows that the perturbative method was a good 
approximation in this case. 

5 Discussions 

We proposed a numerical method to solve the scattering problem of the Dirac equation 
in the background of any CP- violating bubble wall. By use of it, one can evaluate the 
reflection coefficients of chiral fermions to any precision. This will be a powerful tool to 
analyze electroweak baryogenesis based on the charge transport scenario. 

Whether the charge transport or the spontaneous mechanism is effective will be de- 
termined by the dynamics of the EWPT. Recent study of the standard model at finite 
temperature suggests that the first order nature of the EWPT is stronger than that 
expected from the one-loop calculation of the effective potential. Even though efficient 
electroweak baryogenesis requires some extension of the Higgs sector to incorporate spa- 
tially varying CP phase, this fact would mean that the bubble wall would be thinner than 
naively expected from the lowest order effective potential, so that the charge transport 
scenario would play a crucial role in baryogenesis. 

The profiles of the bubble wall which we took in this paper have the kink- type p(x). 
In fact, they should be determined by the dynamics, and we have shown that the modulus 
of the Higgs could be distorted from the kink, when 9(x) is not small for the choice of 
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the effective potential which has a kink-type p(x) for 9 = 0J7|]. For such a profile, we can 
evaluate the flux precisely, by taking sufficiently small a. 
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Figure Captions 

Fig.l: Contour plot of the chiral charge flux, log 10 [— Fq/(uT 3 (Q l — Qr))] for the example 
(i) with O = -1- Here we take u = 0.1 and T = lOOGeV. 

Fig.2: Contour plot of the chiral charge flux, normalized as log 10 [— Fq/(uT 3 (Ql — Qr))] 
for the example (ii), given in Fig. 2 of 0. Here we take the same u and T as Fig. 1. 

Fig.3: AR for the profile given in Fig. 3 of ||. Here we take £ = 0.4 and m = 1. The 
solid line represents the numerical result obtained here. The dashed line stands for 
the previous result based on the perturbative method. 
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